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Abstract 

Using a recently developed thermodynamic integration method, we compute the precise values 
of the excess Gibbs free energy (G e ) of the high density liquid (HDL) phase with respect to the 
crystalline phase at different temperatures (T) in the supercooled region of the Stillinger- Weber 
(SW) silicon [F. H. Stillinger and T. A. Weber, Phys. Rev. B. 32, 5262 (1985)]. Based on the slope 
of G e with respect to T, we find that the absolute entropy of the high density liquid (HDL) phase 
shows a nonmonotonic dependence on temperature at the liquid-liquid transition temperature of 
Tll = 1060 K. Our result is consistent with the earlier observation of a nonmonotonic dependence 
of the enthalpy on temperature in molecular dynamics simulations starting in the HDL phase at a 
temperature just above Tll [S. Sastry and C. A. Angell, Nat. Mater. 2, 739 (2003)]. Our result 
provides a thermodynamic route by means of which the liquid-amorphous transition occurs in SW 
silicon, and possibly, in real silicon. 

PACS numbers: 64.70.Ja, 65.40.gd, 64.70.D-, 05.70.-a 
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I. INTRODUCTION 



The liquid- amorphous transition in silicon, modeled by the Stillinger- Weber (SW) po- 
tential, 1 has been intensely studied^- with an aim of understanding the phase behavior of 
real silicon. In the initial molecular dynamics (MD) studies on SW silicon,-^ it was found 
that the high density liquid (HDL) phase, at a sufficiently slow cooling rate, undergoes a 
sudden transition to a low density amorphous phase at around 1060 K. The nature of the 
low density phase (i.e., whether it is a solid or a liquid) below the transition temperature 
was however not clear. In 2003, Sastry and Angell, through precise and careful measure- 
ments of the diffusivity in MD simulations, unequivocally showed that a low density liquid 
(LDL) phase exists below Tll = 1060 K and hence the transition should be characterized 
as a liquid-liquid transition.- It was also demonstrated that in constant pressure-constant 
enthalpy (NPH) MD simulations starting from the HDL phase at T > 1060 K, the enthalpy 
shows a nonmonotonic dependence on temperature, ultimately leading to the formation of 
the LDL phase.- This was attributed to the release of latent heat during the phase trans- 
formation from the HDL phase to the LDL phase.- The purpose of this work is to perform 
a thermodynamic analysis pertaining to this interesting observation and reconcile all the 
available facts including the earlier observations of a sudden transition at 1060 K in MD 
cooling experiments and the uni-directionality of the transition with respect to the temper- 
ature. This last aspect refers to the finding that when the LDL phase starting at T < Tll 
is heated slowly in MD simulations, it does not convert back to the HDL phase at T > Tll, 
but undergoes crystallization.-^ 

The primary tool that we employ is a recently developed thermodynamic integration 
method^— which enables the precise measurement of the excess Gibbs free energy (G e ) of 
the HDL phase with respect to the crystalline phase. The slope of G e with respect to the 
temperature T yields the excess entropy S e of the HDL phase. We then compare the changes 
in S e with the changes in the absolute entropy of the crystal phase as the temperature is 
lowered to Tll and conclude that the absolute entropy of the HDL phase SjjDL sh° ws a 
nonmonotonic behavior with respect to the temperature at Tll- In what follows, we describe 
the details of our computational method. 
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II. SIMULATIONS OF THE HDL PHASE 



To compute the excess free energy (to be described in the next section), it is important to 
correctly determine the average properties of the HDL phase in the supercooled region. To 
this end, we studied the properties of the supercooled HDL phase by performing isothermal- 
isobaric (NPT) Monte Carlo (MC) simulations at and above Tll = 1060K and zero pressure. 
All of our computations were performed with a system of N = 512 particles in a cubic 
simulation box under periodic boundary conditions. When the free energy barrier preventing 
crystallization of the metastable HDL phase is sufficiently high (T > 1086 K in the present 
case), there is no ambiguity in the computation of the average properties, i.e., all MC 
trajectories starting from different initial configurations converge to the free energy minimum 
relatively quickly. However when the free energy minimum corresponding to the metastable 
HDL phase is shallow, the probability of crossing the free energy barrier is higher and 
therefore proper equilibration is important to compute the correct average properties. Thus, 
in the temperature range of 1060-1082 K, we generated several independent MC trajectories 
by starting from different initial configurations. The equilibration of the trajectory in the 
HDL phase region was ensured by absence of a systematic decrease in the potential energy. 
The equilibrated trajectories at 1060 K and 1075 K are shown in Figs.[I]and[2j respectively. A 
systematic decrease in the potential energy per particle </>', as seen in these figures, indicates 
that the system has crossed the free energy barrier. At 1060 K, the crossing of the free energy 
barrier can be seen more clearly by the systematic decrease in the cumulative average of 0' as 
seen in the inset of Fig. [TJ The region of the trajectory corresponding to the HDL phase can 
thus be identified. The length of the trajectory corresponding to the HDL phase is 68 million 
MC steps at 1075 K (Fig. [2]) and 11 million MC steps at 1060 K (Fig. d]) (Each MC step, 
on average, consisted of two volume change moves and N particle displacement attempts). 
This indicates that the free energy barrier is reduced considerably as the temperature is 
decreased from 1075 K to 1060 K. Further, we find that equilibrated trajectories starting 
from different initial configurations at a given temperature yield the same average properties, 
thus indicating the existence of a local equilibrium corresponding to the metastable HDL 
phase. 

The computed average properties of the HDL phase are reported in Table 1. (The quan- 
tities reported throughout this work are expressed in units of SW potential parameters^ a 
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and e, unless otherwise noted explicitly). The average densities (see Table 1) of our HDL 
phases at 1060 K and 1082 K are in good agreement with the previously reported values 13 
of 0.479 and 0.482 a -3 , respectively. The computed average densities at other temperatures 
agree well with those found in the MD cooling experiments of Beaucage and Mousseau-. 

At 1060 K, the MC trajectory upon crossing the free energy barrier, stabilizes in the 
LDL phase (see Fig. [Tj), while at higher temperatures (see Fig. [2]), the MC trajectory yields 
a crystalline phase with defects (d-crystal phase). The LDL phase at 1060 K (see Fig. [1]) has 
a average total energy (~ — 1.793e per particle or 375 kJ/mol) comparable to the enthalpy 
of the LDL phase at 1055 K found in Ref. O (see Fig. 1 of Ref. 6). The transition from 
the HDL phase to the LDL or the d-crystal phases at a given temperature involves entirely 
homogeneous states. This is deduced based on the observation that the following fluctuation 
relation 14 is satisfied by the MC trajectory, 

({p v - (p v ))(v - <y))) = -t, (l) 

where P v is the instantaneous pressure in the NPT-MC simulations calculated from the virial 
relation^ at the given instantaneous volume V and T is the externally set temperature. The 
symbol (■ ■ •) represents the average taken over the entire trajectory of the isothermal isobaric 
MC simulations. This fluctuation relation is derived by assuming^ 4 - that the instantaneous 
fluctuations represents a change in state from one homogeneous phase to the other. The 
relation is satisfied by the entire trajectory consisting of the HDL, LDL and the defect crystal 
regions, as well as by the partial trajectories consisting only of the HDL phase region. 
Formally, the partition function of a metastable state is defined as follows: 16 

Y(T, P, N) = J dV J dp N J dr N exp[-BH - f3PV - pw{r x , • • • , r N )\, (2) 

where H is the system Hamiltonian, P is the external pressure, and w = w(r x , ■ ■ ■ , tn) 
is an additional potential energy (to be imposed externally) that reflects back the system 
from the top of the free energy barrier towards the free energy minimum.— It is zero in the 
single phase region and takes arbitrarily large and positive values in regions beyond the free 
energy barrier. Thus the role of the constraining potential w is to confine the system to the 
metastable free energy minimum.— The probability density in isothermal isobaric ensemble 
is then given by, 

/ at ^ T „ exp(-/W - BH - Bw) 



Applying the Gibbs entropy formula S = — ks(logp), we have 



(4) 



Since the Gibbs free energy is defined as G 
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Comparing Eqs. (@| and (jSJ), we get 
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When the free energy barrier is large, the probability of reaching the top of the barrier is 
small and hence the quantity (w) in Eq. flj]) can be neglected. However, when the free energy 
barrier is small, the ensemble average (w) cannot be neglected. In such a case, based on 
Eq. fllj) we can conclude that G ^ (H)—TS+P(V). In all of our NPT-MC simulations, we do 
not impose the constraint function w [see Eq. fl2}], but we can still obtain the correct average 
properties of the HDL phase (e.g., (<j)') and (p)) by ensuring equilibration of the trajectory. 
This is feasible due to the presence of a free energy barrier at T > 1060 K. Nonetheless, the 
consideration of the reflective barrier w is, in principle, essential for defining the partition 
function of the metastable HDL phase.— 

III. COMPUTATION OF EXCESS GIBBS FREE ENERGY 

We computed the excess Gibbs free energy difference G e = Gjjql — ^crystal between 
the HDL and the crystal phases, by applying the constrained fluid A integration method- 
in the isothermal-isobaric ensemble.— Recently, the method was found to predict the 
melting point of SW silicon accurately— This is a thermodynamic integration method in 
which the liquid and the crystal phase are connected directly through a 3-stage reversible 
path. In stage 1 of the reversible path, which starts from the liquid phase, the strength of 
the interaction potential is reduced linearly so that the system approaches an ideal gas-like 
behavior. The expression for the potential energy in this stage is given by,- 1 ^ 



where r] is a constant that determines the effective strength of the interaction potential at 
the end of stage 1 and U is the original inter-particle potential ( SW potential, in the present 



1 (A 1 ) = (l-r / A 1 )[/, 



(7) 
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case). The parameter Ai defines the states along the path and varies from to 1. As the 
system becomes less attractive it tends to expand. However, to maintain the reversibility of 
the path in stages 2 and 3, it is necessary to impose a maximum volume constraint. 10 The 
maximum constrained volume V m is chosen such that it is slightly larger than the average 
volumes of the liquid and the solid phase (whichever is larger). At the same time, V m should 
not affect the free energies of either of these phases.— Such a volume can be straightforwardly 
chosen based on the histogram of volume fluctuations for the two phases. As in an earlier 
work on SW silicon,— we chose V m to correspond to a density of 0.4<t~ 3 , i.e., V m = N/OA. 
At the end of stage 1, we get a compressed gas phase due to the constraint on the maximum 
volume.— This process is depicted pictorially in Fig. [3J 

In stage 2, we force the particles to form a crystalline structure by imposing an external 
potential in the form of Gaussian potential wells distributed on the ideal crystal lattice.- 
The strength of the inter-particle potential energy is held fixed during this stage. The total 
potential energy for stage 2 is given by^ 1 ^ 



The Gaussian external potential imposed during this and the subsequent stage is given by 
U ex t = J2i J2k aex P(~brf k )- 9 Here, the index 'i' refers to the system particle and the index 'k' 
is a Gaussian potential well. The Gaussian well does not act on a specific particle, but exerts 
an influence over all the particles in its vicinity.- 1 ^ The values of the Gaussian parameters 
are taken to be the same as in the earlier work:— r\ = 0.9, a = — 1.892e and b = 8.0a. These 
values are so chosen that the constrained crystalline state obtained at the end of stage 2 has 
almost the same energy and density as the desired crystal phase.— ^ 

In stage 3, the Gaussian external potential is reduced linearly to zero, while the strength 
of the potential energy is increased linearly to its original value.- The potential energy 
expression for this stage is given by^ 



At the end of this stage, we get the desired crystalline phase as shown pictorially in Fig. [3j 

th 

The Gibbs free energy change for the i stage of the path can be obtained by numerical 
integration : 



2 (A 2 ) = (i - v )u + \ 2 u ext . 



(8) 



<f> 3 (X 3 ) = [(l-r l ) + \ S 7j\U+(l-\ 3 )U t 



ext • 



(9) 




(10) 
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where (• • •) represents the isothermal-isobaric ensemble average at a given value of Aj. The 
integrands for the various stages of the path at 1065 K and zero pressure are plotted in 
Figs. HHSJ It can be seen from these figures that the value of the integrand for the forward 
and the reverse paths agree well at each point. This shows that there is no hysteresis present 
along the path, as found earlier.— 

Throughout this work, we have used the Bennett Acceptance Ratio (BAR) method,— to 
compute the Gibbs free energy between the adjacent states along the entire path. Accord- 
ing to the BAR method,— the Gibbs free energy difference AG = G% — Go between two 
equilibrium states '0' and '1', for a given value of the constant C, is given by the following 
equation:— 

— — = log — — + G-log— , 11 

k B T Eo /Wi - P0o - C) n 

where f(x) = 1/(1 + e x ) is the Fermi function, an d J2i represent the sums over Fermi 
functions sampled in '0' and '1' ensembles, respectively. The total potential energies in the 
two ensembles are represented by 0o and 0i in Eq. (iTTj) and the total number of samples of 
the perturbation energies (or equivalently the Fermi functions) collected in MC simulations 
in the two ensembles are no and n\. In principle, the above equation yields the correct 
value of AG for any value of the constant G. In practice, due to the limited computational 
power, we cannot sample the perturbation energies in the two ensembles over all possible 
configurations. Thus, Bennett showed that the optimum value of G, which yields minimum 
error in the estimation of AG, is given by^ 

—— = G-log— . 12 

k B T n 

Combining the above two equations, we gelr^ 

J2 Wo - Pfa + C) = £ /(/30i - /30o - C). (13) 
1 o 

The value of the G satisfying the above equation is substituted in Eq. (fl2|) to yield the 



optimal estimate of AG. As mentioned in Ref. [18|, the accuracy of the above method 
depends on the degree of overlap in the configuration space between the two ensembles. The 
larger the value of the two sums appearing in Eq. ( lT3"j) . the greater is the configuration space 
overlap. As prescribed by Bennett, the sum values should be much greater than unity.— 

In all our computations, we ensured that the sum values are of the order of 10 5 -10 6 , by (i) 
performing simulations at Aj values that are sufficiently close to each other (see Figs.HHB]) and 



(ii) performing sufficiently long simulation runs at each of the Aj values. The perturbation 
energies [(</>i — <fio) or (0 O — <j>i) in Eq. (|13|) ] required in the BAR method were collected after 
every MC step. At Ai = 0.0 in stage 1, we used the properly equilibrated trajectories of the 
HDL phases, as described in the Sec. II, to collect the perturbation energy data. In stage 
1, we performed upto 20 million MC steps from Ai = to 0.4. In the region from Ai = 0.4 
(stage 1) to A3 = 0.99 (stage 3), we performed 0.4-0.8 million MC steps. The main source of 
statistical error was found to be towards the end of stage 3 due to the center of mass motion 



of the crystal phase as explained in detail in Ref. [12j. To address this problem, we performed 
up to 10 million MC steps from A3 = 0.992 to A3 = 0.999 to collect the perturbation energy 
data. In the last interval of stage 3 from A3 = 0.999 to 1, there is a large change in the 
integrand (due to center of mass motion of the crystalline phase) as seen in the inset of 
Fig. [61 In order to minimize the statistical error, upto 30 million MC steps were performed 
at A 3 = 1. We ensured that the value of the sums J2o an d Xa [ see Eq. (|T3|) ] over the Fermi 
function is above 10 5 for all the intervals towards the end of stage 3. Further, in order to 
improve the accuracy three to four independent simulation runs along the entire path at all 

the temperatures. The statistical error for the ith stage was computed by using the formula 

1 11 

Oi = a(AGj)/n/ , where cr(AGj) is the standard deviation in the AGi value and rii is the 
number of statistically independent measurements. The total statistical error was computed 
by adding the errors for the individual stages. 

IV. EXCESS AND ABSOLUTE ENTROPY OF THE HDL PHASE 

Figure 7 shows the values of negative excess Gibbs free energy — G e computed at various 
temperatures (also see Table 1) in the supercooled region at and above Tn. The G e values 
are comparable to those calculated by Broughton and L>^ (~ 25. 3e for N = 512 at 1060 K 



obtained by linear interpolation of the chemical potential data in Table III of Ref. (2). The 
excess entropy of the HDL phase S e is equal to the slope of the curve, which we compute 
by the forward difference method: 

s e = - * +1 y , (14) 

where T i+ i > T{. The values of S e computed by above equation (for N = 512) are reported 
in Table 1. The forward difference method gives the average slope (i.e. average value of S e ) 
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in the temperature interval (T,T + i). From Table 1, we find that the average value of S e in 
the interval (1060 K,1065 K) is higher by at least 250k B (after considering the error bars) 
compared to that in the interval (1065 K, 1070 K). Thus, as the temperature is decreased 
by AT = —5 K from 1067.5 to 1062.5 K (the mid-points of above two intervals) the change 
in the excess entropy is AS* 6 > 250 k B - 

To estimate the change in entropy of the crystalline phase, we evaluated its constant 
volume heat capacity at 1060 K by using the following relation 19 

Cv 3 {{84?) 

Nk B 2 N(k B T) 2 ' { ' 

where the quantity 5(f) = <fi — {(f)) is the instantaneous fluctuation in the total potential 
energy of the system. The average over the potential energy fluctuations was computed in 
the constant volume-constant temperature (NVT) MC simulations of the crystal phase and 
substituting this value in the above equation yields Cy = 3.7 Nk B at 1060 K. We neglect 
the difference between the heat capacities of the crystal phase at constant pressure and at 
constant volume, i.e., Cp ~ Cy = 3.7Nk B . Then the change in absolute entropy of the 
crystal phase can be estimated by the following formula, 

(T + AT\ 

A£«C P ]og(-^— J, (16) 

where we assume Cp to be constant over the temperature range of AT. Using the above 
equation, we find that the change in the entropy A^j-yg^^ of the crystal phase is approx- 
imately —9k B (for N = 512) when the temperature is reduced from 1067.5 K to 1062.5 K 
(AT = —5 K). Since AS e = ASjjjjl ~~ ^^crystab we nnc ^ that ^^HDL — 259 k B for the 
change in temperature AT = —5 K from 1067.5 K to 1062.5 K. In other words, the HDL 
phase entropy shows a nonmonotonic dependence on temperature near 1060 K. 

This is the main result of our work and we discuss its implications in the next section. 



|P) <0. (17) 



V. SUMMARY AND CONCLUSIONS 



In this work, we have computed precisely the excess Gibbs free energy G e of the HDL 
phase at different temperatures by using a recently developed thermodynamic integration 
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method.- 1 ^ The slope of the G e versus T curve yields the excess entropy of the HDL phase 
with respect to the crystal phase. Based on the change in the slope of the curve with respect 
to the temperature, we conclude that the absolute entropy of the HDL phase increases as 
the temperature is decreased to 1060 K. For a stable equilibrium phase, the entropy cannot 
increase as the temperature decreases. However, our result applies to the metastable HDL 
phase with a shallow free energy minimum. In NPT MC simulations starting from the 
HDL phase, we observe (see Sec. II) that system undergoes crystallization after a certain 
number (e.g., 11 million at 1060 K) of MC steps. Therefore the state of the isolated system 
containing the HDL phase and the heat bath must be considered as a constrained equilibrium 
state. The constraint may be imposed by means of the additional potential energy function 
w(ri, ■ ■ • , r N ) as given in Eq. fl2J).— 

Our findings are consistent with the nonmonotonic dependence of enthalpy on tempera- 
ture observed by Sastry and Angell. In the NPH MD simulations,- the nonmonotonic loop 
starts in the HDL phase at a temperature just above 1060 K. As the enthalpy (or equiv- 
alently the internal energy) is reduced at zero pressure, no mechanical work is performed, 
and hence the lowering of internal energy must be due to the removal of heat, which im- 
plies reduction of the entropy. According to our computations, the absolute entropy of the 
HDL phase increases as the temperature is reduced to 1060 K. Therefore, the HDL phase at 
T > 1060 K, upon decreasing the enthalpy, cannot transform into the HDL phase at a lower 
temperature. Thus the only alternative is to transform to a phase (which is intermediate 
between the HDL and the LDL) at a temperature equal to or greater than the starting 
HDL phase temperature. Note that the LDL phase has a lower enthalpy (see Fig. [T]) and is 
expected to have a lower entropy than the HDL phase since it is structurally closer to the 
crystalline phase.-^i Thus, our computations show that the nonmonotonic loop starting in 
the HDL phase at T > 1060 K, 6 is induced by the nonmonotonic dependence of the HDL 
phase entropy on temperature. 

In the previous MD studies^^ it was observed that the HDL phase, at a sufficiently slow 
cooling rate, transforms into a low density amorphous phase at or below 1060 K. (Here by 
amorphous phase, we do not necessarily mean the LDL phase, but any intermediate phase 
between the HDL and the LDL phases). Since the process of cooling at zero pressure nec- 
essarily involves reduction of the entropy, the HDL phase at T > 1060 K will not transform 
into HDL phase at 1060 K, since the latter has a higher entropy. Instead, the transforma- 
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tion to a low density amorphous phase (with a lower entropy) at 1060 K or below is the 
most natural outcome. Thus, the liquid-amorphous transition observed in the MD cooling 
experiments 4 - 1 ^ can be explained on the basis of our result. 

Recent first principles MD simulations of supercooled silicon^ have found a liquid-liquid 
(LL) transition of a similar nature to that in the case of SW silicon. Thus we expect that 
SW silicon shows a qualitatively similar behavior to that of the real silicon and hence the 
nonmonotonic dependence of the HDL phase entropy on temperature could be the underlying 
cause of the LL transition in real silicon. 

Finally, if one agrees that the nonmonotonic dependence of the entropy S^DL on tem- 
perature T is the cause of the HDL to LDL (or amorphous) transition as argued above, 
then the reverse transition from the LDL to the HDL phase is feasible only if the LDL 
phase entropy decreases as the temperature is increased above Tll, be., if Sldl a ^ so snows 
a nonmonotonic dependence on T. Since the LDL phase is structurally much closer to the 
stable crystal phase&^i and also has a heat capacity comparable to that of the crystal,— 
its thermodynamic phase behavior should also be similar to the crystal phase and hence 
the above scenario does not seem likely. Hence the uni-directionality of the transition as 
observed earlier-^ can be rationalized on the basis of the expectation that Sldl shows a 
monotonic dependence on T similar to the crystalline phase entropy. 
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TABLE I: The excess Gibbs free energy of the HDL phase with respect to the crystalline phase 
is listed at various temperatures for N = 512 and zero pressure. The excess entropy listed in the 
third column is computed by forward difference method as explained in the text. The G e data 
contained in this table is also plotted in Fig. 7. The fourth and the fifth columns contain the 
average densities and average potential energy per particle of the HDL phases. The statistical 
error in the computation of (ft) is 5 x 1(T 4 e for T < 1082 K, and 1 x 1(T 3 e for T > 1082 K. The 
statistical error in the computation of (p) is 1 x 10~ 3 o --3 at all temperatures listed. 



T(K) 


G e (e) 


S e /k B 


(P) (v~ 3 ) 


W) (e) 


1060 


26.09 ± 0.02 


1257 ± 200 


0.478 


-1.8218 


1065 


25.84 ±0.02 


604 ± 200 


0.478 


-1.8216 


1070 


25.72 ±0.02 


855 ± 200 


0.479 


-1.8195 


1075 


25.55 ±0.02 


683 ± 144 


0.479 


-1.8184 


1082 


25.36 ± 0.02 


1006 ± 252 


0.482 


-1.815 


1086 


25.20 ±0.02 


503 ± 252 


0.482 


-1.814 


1090 


25.12 ±0.02 


930 ± 100 


0.482 


-1.813 


1100 


24.75 ± 0.02 




0.483 


-1.810 
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x 10 6 MC steps 



FIG. 1: The MC trajectory in terms of the potential energy per particle ((f)' = <fi/N) at T = 1060 
K. The squares represent block averages taken after every 0.2 million MC steps. The LDL phase 
region starts (roughly) from the 52 millionth MC step onwards, as shown. The inset shows the 
region of the trajectory corresponding to the HDL phase. The circles in the inset are the cumulative 
averages obtained after the given number of MC steps on the abscissa. The dashed line in the inset 
is the average energy (which is the same as the cumulative average energy) of the HDL phase 
computed from the initial 11 million MC steps of the trajectory. A systematic and continuous 
decrease in the cumulative average signals the crossing of the free energy barrier. 



14 




FIG. 2: The MC trajectory in terms of block averages of potential energy per particle at 1075 
K. The squares are the block averages taken after every 0.2 million MC steps. The HDL phase 
region corresponds to the initial 68 million MC steps of the trajectory, as shown. The dashed 
horizontal line represents the average energy of the HDL phase. The trajectory ultimately leads to 
the formation of the crystalline phase with defects (d-crystal phase). The portion of the trajectory 
corresponding to the d-crystal phase starts (roughly) from 80 millionth MC step onwards. 
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FIG. 3: Schematic diagram of the thermodynamic integration path connecting the liquid and the 
solid states at constant temperature and constant external pressure.— 1 ^ The stops in the piston- 
cylinder arrangement represent the maximum volume constraint as described in the text. 
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FIG. 4: The integrand (in units of the SW potential parameter e) in Eq. (|10p as a function of Ai 
for stage 1 at 1065 K, P = and N = 512. 
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FIG. 5: The same as in Fig. [H but for stage 2. 
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FIG. 6: The same as in Fig.UJ but for stage 3. The inset shows the region of the plot from A3 = 0.98 
to 0.999. 
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FIG. 7: The negative excess Gibbs free energy G e = Gjjdl — C cr y sta j (in units of e) as a function 
of temperature (first and second column of Table 1) for P = and N = 512. The error bars (see 
Table 1) are smaller than the size of the symbols and hence cannot be seen. The slope of the 
curve yields the excess entropy S e of the HDL phase with respect to the crystalline phase. A large 
increase in the slope is seen as the temperature is lowered to 1060 K, indicating an increase in the 
excess and the absolute entropy of the HDL phase as discussed in the text. 
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